HS project

Zoom in IBL and mORS cells

This file is used to analyse the immune cells dataset.

library(dplyr)
library(patchwork)
library(ggplot2)
library(ComplexHeatmap)
library(org.Mm.eg.db)

.libPaths()

## [1] "/usr/local/lib/R/library"

Preparation

In this section, we set the global settings of the analysis. We will store data there :

save_name = "iblmors"
out_dir = "."

We load the dataset :

sobj = readRDS(paste0(out_dir, "/", save_name, "_sobj.rds"))
sobj

## An object of class Seurat 
## 16701 features across 3532 samples within 1 assay 
## Active assay: RNA (16701 features, 2000 variable features)
##  6 dimensional reductions calculated: RNA_pca, RNA_pca_20_tsne, RNA_pca_20_umap, harmony, harmony_20_umap, harmony_20_tsne

We load the sample information :

sample_info = readRDS(paste0(out_dir, "/../../1_metadata/hs_hd_sample_info.rds"))
project_names_oi = sample_info$project_name

graphics::pie(rep(1, nrow(sample_info)),
              col = sample_info$color,
              labels = sample_info$project_name)

Here are custom colors for each cell type :

color_markers = readRDS(paste0(out_dir, "/../../1_metadata/hs_hd_color_markers.rds"))

data.frame(cell_type = names(color_markers),
           color = unlist(color_markers)) %>%
  ggplot2::ggplot(., aes(x = cell_type, y = 0, fill = cell_type)) +
  ggplot2::geom_point(pch = 21, size = 5) +
  ggplot2::scale_fill_manual(values = unlist(color_markers), breaks = names(color_markers)) +
  ggplot2::theme_classic() +
  ggplot2::theme(legend.position = "none",
                 axis.line = element_blank(),
                 axis.title = element_blank(),
                 axis.ticks = element_blank(),
                 axis.text.y = element_blank(),
                 axis.text.x = element_text(angle = 30, hjust = 1))

This is the projection of interest :

name2D = "harmony_20_tsne"

Visualization

Gene expression

We visualize gene expression for some markers :

features = c("percent.mt", "percent.rb", "nFeature_RNA")

plot_list = lapply(features, FUN = function(one_gene) {
  Seurat::FeaturePlot(sobj, features = one_gene,
                      reduction = name2D) +
    ggplot2::theme(aspect.ratio = 1) +
    ggplot2::scale_color_gradientn(colors = aquarius::color_gene) +
    Seurat::NoAxes()
})

patchwork::wrap_plots(plot_list, ncol = 3)

Clusters

We visualize clusters :

cluster_plot = Seurat::DimPlot(sobj, reduction = name2D, label = TRUE) +
  Seurat::NoAxes() +
  ggplot2::theme(aspect.ratio = 1)
cluster_plot

Cell type

We visualize cell type split by sample :

plot_list = aquarius::plot_split_dimred(sobj,
                                        reduction = name2D,
                                        split_by = "project_name",
                                        group_by = "cell_type",
                                        split_color = setNames(sample_info$color,
                                                               nm = sample_info$project_name),
                                        group_color = color_markers,
                                        bg_pt_size = 0.5, main_pt_size = 0.5)

plot_list[[length(plot_list) + 1]] = cluster_plot

patchwork::wrap_plots(plot_list, ncol = 4) &
  Seurat::NoLegend()

Cluster type

We summarize major cell type by cluster :

cell_type_clusters = sobj@meta.data[, c("cell_type", "seurat_clusters")] %>%
  table() %>%
  prop.table(., margin = 2) %>%
  apply(., 2, which.max)
cell_type_clusters = setNames(levels(sobj$cell_type)[cell_type_clusters],
                              nm = names(cell_type_clusters))

We define cluster type :

sobj$cluster_type = cell_type_clusters[sobj$seurat_clusters] %>%
  as.factor()
table(sobj$cluster_type, sobj$cell_type)

##      
##       CD4 T cells CD8 T cells Langerhans cells macrophages B cells cuticle
##   IBL           0           0                0           0       1       1
##   ORS           2           0                0           0       3       7
##      
##       cortex medulla  IRS proliferative  IBL  ORS  IFE HFSC sebocytes
##   IBL      0       0    0             9 1465    5    7   19         3
##   ORS      2      20    7            37   26 1785   65   48        20

We subset color_markers :

color_markers = color_markers[levels(sobj$cluster_type)]

We compare cluster annotation and cell type annotation :

p1 = Seurat::DimPlot(sobj, group.by = "cell_type",
                     reduction = name2D, cols = color_markers) +
  ggplot2::labs(title = "Cell type") +
  Seurat::NoAxes() +
  ggplot2::theme(aspect.ratio = 1,
                 plot.title = element_text(hjust = 0.5))

p2 = Seurat::DimPlot(sobj, group.by = "cluster_type",
                     reduction = name2D, cols = color_markers) +
  ggplot2::labs(title = "Cluster type") +
  Seurat::NoAxes() +
  ggplot2::theme(aspect.ratio = 1,
                 plot.title = element_text(hjust = 0.5))

patchwork::wrap_plots(p1, p2, guides = "collect")

Barplot

We make a barplot to compare IBL and ORS populations in HS vs HD samples. The proportion of IBL is indicated on top of bars.

quantif = table(sobj$sample_identifier,
                sobj$cluster_type) %>%
  as.data.frame.table() %>%
  `colnames<-`(c("Sample", "cell_type", "nb_cells")) %>%
  dplyr::group_by(Sample) %>%
  dplyr::mutate(total_cells = sum(nb_cells)) %>%
  as.data.frame() %>%
  dplyr::filter(cell_type == "IBL") %>%
  dplyr::mutate(prop_ibl = nb_cells / total_cells) %>%
  dplyr::mutate(prop_ibl = 100*round(prop_ibl, 4))

sobj$seurat_clusters = factor(sobj$seurat_clusters,
                              levels = names(sort(cell_type_clusters)))

aquarius::plot_barplot(df = table(sobj$sample_identifier,
                                  sobj$seurat_clusters) %>%
                         as.data.frame.table() %>%
                         `colnames<-`(c("sample_identifier", "clusters", "nb_cells")),
                       x = "sample_identifier", y = "nb_cells", fill = "clusters",
                       position = position_fill()) +
  ggplot2::scale_fill_manual(values = c(colorRampPalette(c("chartreuse1", "chartreuse4"))(table(sort(cell_type_clusters))["IBL"]),
                                        colorRampPalette(c("cadetblue1", "cadetblue4"))(table(sort(cell_type_clusters))["ORS"])),
                             breaks = names(sort(cell_type_clusters)),
                             name = "Cell type") +
  ggplot2::geom_label(data = quantif, inherit.aes = FALSE,
                      aes(x = .data$Sample, y = 1.05, label = .data$prop_ibl),
                      label.size = 0, size = 5)

Differential expression

In this section, we perform DE between :

• inner bulge layer (IBL) and outer root sheath (ORS) populations
• healthy donors (HD) and HS patients (HS) within ORS population
• healthy donors (HD) and HS patients (HS) within IBL population
• cluster 5 vs rest of ORS, because this cluster is specific to HS compared to HD

We save the results in a list :

list_results = list()

IBL vs ORS

We change cell identities to cluster type :

Seurat::Idents(sobj) = sobj$cluster_type

table(Seurat::Idents(sobj), sobj$sample_type)

##      
##         HS   HD
##   IBL 1284  226
##   ORS 1563  459

We identify specific markers for each population :

mark = Seurat::FindMarkers(sobj, ident.1 = "IBL", ident.2 = "ORS")

mark = mark %>%
  dplyr::filter(p_val_adj < 0.05) %>%
  dplyr::arrange(-avg_logFC, pct.1 - pct.2)

list_results[["IBL_vs_ORS"]] = mark

dim(mark)

## [1] 1027    5

head(mark, n = 20)

##                  p_val avg_logFC pct.1 pct.2     p_val_adj
## KRT16     0.000000e+00  4.323925 0.920 0.116  0.000000e+00
## KRT6B     0.000000e+00  3.627468 0.921 0.285  0.000000e+00
## FABP5     0.000000e+00  3.348761 0.993 0.509  0.000000e+00
## KRT17     0.000000e+00  3.020736 0.987 0.733  0.000000e+00
## KRT6A    1.267920e-155  2.994730 0.764 0.476 2.117554e-151
## CST6     1.834243e-160  2.741295 0.391 0.038 3.063369e-156
## KRT6C    9.379153e-283  2.653075 0.568 0.038 1.566412e-278
## TMSB4X    0.000000e+00  2.356496 0.989 0.864  0.000000e+00
## LGALS1    0.000000e+00  2.320083 0.929 0.484  0.000000e+00
## S100A2    0.000000e+00  2.134213 0.996 0.977  0.000000e+00
## CALML3    0.000000e+00  2.003586 0.972 0.686  0.000000e+00
## GJA1      0.000000e+00  1.973282 0.875 0.586  0.000000e+00
## GJB6      0.000000e+00  1.907928 0.916 0.656  0.000000e+00
## SBSN     9.451796e-242  1.751399 0.534 0.047 1.578544e-237
## APOE      0.000000e+00  1.744652 0.983 0.772  0.000000e+00
## NDUFA4L2 1.058296e-239  1.730465 0.697 0.217 1.767461e-235
## LYPD3     0.000000e+00  1.680352 0.738 0.135  0.000000e+00
## TUBB2A    0.000000e+00  1.667800 0.756 0.181  0.000000e+00
## SDC1      0.000000e+00  1.662321 0.838 0.517  0.000000e+00
## GJB2     1.288894e-254  1.634612 0.826 0.595 2.152582e-250

There are 1027 genes differentially expressed. How many for each population ?

# IBL
mark %>%
  dplyr::filter(avg_logFC > 0) %>%
  nrow()

## [1] 434

# ORS
mark %>%
  dplyr::filter(avg_logFC < 0) %>%
  nrow()

## [1] 593

We represent the information on a figure :

mark$gene_name = rownames(mark)
mark_to_label = rbind(
  # up-regulated in IBL
  mark %>% dplyr::top_n(., n = 20, wt = avg_logFC),
  # up-regulated in ORS
  mark %>% dplyr::top_n(., n = 20, wt = -avg_logFC),
  # representative and selective for IBL
  mark %>% dplyr::top_n(., n = 20, wt = (pct.1 - pct.2)),
  # representative and selective for ORS
  mark %>% dplyr::top_n(., n = 20, wt = -(pct.1 - pct.2))) %>%
  dplyr::distinct()

ggplot2::ggplot(mark, aes(x = pct.1, y = pct.2, col = avg_logFC)) +
  ggplot2::geom_abline(slope = 1, intercept = 0, lty = 2) +
  ggplot2::geom_point() +
  ggrepel::geom_label_repel(data = mark_to_label,
                            aes(x = pct.1, y = pct.2, label = gene_name),
                            col = "black", fill = NA, size = 3, label.size = NA) +
  ggplot2::labs(title = "Differentially expressed genes",
                subtitle = "between IBL and ORS") +
  ggplot2::scale_color_gradient2(low = aquarius::color_cnv[1],
                                 mid = aquarius::color_cnv[2],
                                 high = aquarius::color_cnv[3],
                                 midpoint = 0) +
  ggplot2::theme_classic() +
  ggplot2::theme(plot.title = element_text(hjust = 0.5),
                 plot.subtitle = element_text(hjust = 0.5))

We represent the best genes on a violin plot, group by cluster type and split by sample type :

mark_to_label = mark_to_label %>%
  dplyr::arrange(pct.1 - pct.2)

plot_list = lapply(mark_to_label$gene_name, FUN = function(gene) {
  p = Seurat::VlnPlot(sobj, features = gene, pt.size = 0.001,
                      group.by = "cluster_type",
                      split.by = "sample_type") +
    ggplot2::scale_fill_manual(breaks = c("HS", "HD"),
                               values = c("#C55F40", "#2C78E6")) +
    ggplot2::theme(axis.title.x = element_blank(),
                   axis.title.y = element_blank(),
                   legend.position = "none")
  return(p)
})

patchwork::wrap_plots(plot_list, ncol = 5)

On the figure, we can see that some specific markers for a population are indeed specific to a sample type rather than the whole population.

Cluster 5 within ORS

We subset the Seurat object and change cell identities to sample type.

subsobj = subset(sobj, cluster_type == "ORS")
subsobj$seurat_clusters = base::droplevels(subsobj$seurat_clusters)
Seurat::Idents(subsobj) = subsobj$seurat_clusters

table(subsobj$seurat_clusters)

## 
##   0   1   3   5   7 
## 669 438 384 302 229

We identify specific markers for each population :

mark = Seurat::FindMarkers(subsobj, ident.1 = 5)

mark = mark %>%
  dplyr::filter(p_val_adj < 0.05) %>%
  dplyr::arrange(-avg_logFC, pct.1 - pct.2)

list_results[["cluster5_vs_ORS"]] = mark

dim(mark)

## [1] 395   5

head(mark, n = 20)

##                  p_val avg_logFC pct.1 pct.2     p_val_adj
## SLPI     7.619364e-128 1.5306935 0.765 0.164 1.272510e-123
## MT1X      4.989262e-73 1.4231367 0.954 0.680  8.332567e-69
## KRT15    2.275685e-123 1.3705217 1.000 0.874 3.800622e-119
## LY6D      1.984327e-88 1.2250993 0.503 0.080  3.314025e-84
## IGFBP3   1.494292e-140 1.1855371 0.437 0.016 2.495618e-136
## IFI27     1.438175e-99 1.1442636 0.556 0.085  2.401895e-95
## MT1E      2.040603e-70 1.0839955 0.957 0.659  3.408012e-66
## WNT3     1.095829e-212 1.0819468 0.924 0.134 1.830143e-208
## AQP3      5.658795e-83 1.0173107 0.977 0.733  9.450753e-79
## NEAT1     3.110435e-52 0.9735049 0.990 0.905  5.194737e-48
## TSC22D3   1.298633e-75 0.9601088 0.742 0.258  2.168848e-71
## MTRNR2L1  4.400300e-32 0.9329154 0.815 0.593  7.348942e-28
## CXCL14    7.450833e-63 0.9282610 0.997 0.840  1.244364e-58
## ZFP36     1.579077e-42 0.9002595 0.742 0.406  2.637216e-38
## GLUL      7.969325e-58 0.8930699 0.934 0.699  1.330957e-53
## ID1       2.365404e-17 0.8903797 0.652 0.442  3.950461e-13
## FOS       1.096048e-36 0.8789679 0.967 0.837  1.830511e-32
## IL1R2    5.014612e-122 0.8759829 0.675 0.110 8.374904e-118
## IL18     6.622739e-101 0.8511305 0.970 0.456  1.106064e-96
## HOPX      1.939070e-75 0.8395659 0.871 0.417  3.238441e-71

There are 395 genes differentially expressed. How many for each population ?

# Cluster 5
mark %>%
  dplyr::filter(avg_logFC > 0) %>%
  nrow()

## [1] 225

# Other ORS
mark %>%
  dplyr::filter(avg_logFC < 0) %>%
  nrow()

## [1] 170

We represent the information on a figure :

mark$gene_name = rownames(mark)
mark_cluster5 = rbind(
  # up-regulated in cluster 5
  mark %>% dplyr::top_n(., n = 20, wt = avg_logFC),
  # up-regulated in other ORS
  mark %>% dplyr::top_n(., n = 20, wt = -avg_logFC),
  # representative and selective for cluster 5
  mark %>% dplyr::top_n(., n = 20, wt = (pct.1 - pct.2)),
  # representative and selective for other ORS
  mark %>% dplyr::top_n(., n = 20, wt = -(pct.1 - pct.2))) %>%
  dplyr::distinct()

ggplot2::ggplot(mark, aes(x = pct.1, y = pct.2, col = avg_logFC)) +
  ggplot2::geom_abline(slope = 1, intercept = 0, lty = 2) +
  ggplot2::geom_point() +
  ggrepel::geom_label_repel(data = mark_cluster5,
                            aes(x = pct.1, y = pct.2, label = gene_name),
                            col = "black", fill = NA, size = 3, label.size = NA) +
  ggplot2::labs(title = "Differentially expressed genes",
                subtitle = "between cluster 5 and other ORS") +
  ggplot2::scale_color_gradient2(low = aquarius::color_cnv[1],
                                 mid = aquarius::color_cnv[2],
                                 high = aquarius::color_cnv[3],
                                 midpoint = 0) +
  ggplot2::theme_classic() +
  ggplot2::theme(plot.title = element_text(hjust = 0.5),
                 plot.subtitle = element_text(hjust = 0.5))

We represent a subset of genes on a heatmap :

features_oi = mark %>%
  dplyr::filter(p_val_adj < 0.05) %>%
  dplyr::filter(abs(avg_logFC) > 0.5) %>%
  rownames()

length(features_oi)

## [1] 112

We prepare the scaled expression matrix :

mat_expression = Seurat::GetAssayData(subsobj, assay = "RNA", slot = "data")[features_oi, ]
mat_expression = Matrix::t(mat_expression)
mat_expression = dynutils::scale_quantile(mat_expression) # between 0 and 1
mat_expression = Matrix::t(mat_expression)
mat_expression = as.matrix(mat_expression) # not sparse

dim(mat_expression)

## [1]  112 2022

We prepare the heatmap annotation :

ha_top = ComplexHeatmap::HeatmapAnnotation(
  cluster_type = subsobj$cluster_type,
  sample_type = subsobj$sample_type,
  cluster = subsobj$seurat_clusters,
  col = list(cluster_type = color_markers,
             sample_type = setNames(nm = c("HS", "HD"),
                                    c("#C55F40", "#2C78E6")),
             cluster = setNames(nm = levels(subsobj$seurat_clusters),
                                aquarius::gg_color_hue(length(levels(subsobj$seurat_clusters))))))

And the heatmap :

ht = ComplexHeatmap::Heatmap(mat_expression,
                             col = aquarius::color_cnv,
                             # Annotation
                             top_annotation = ha_top,
                             # Grouping
                             column_title = NULL,
                             cluster_rows = TRUE,
                             cluster_columns = TRUE,
                             show_column_names = FALSE,
                             # Visual aspect
                             show_heatmap_legend = TRUE,
                             border = TRUE)

ComplexHeatmap::draw(ht,
                     merge_legend = TRUE,
                     heatmap_legend_side = "bottom",
                     annotation_legend_side = "bottom")

IBL : HS vs HD

We subset the Seurat object and change cell identities to sample type.

subsobj = subset(sobj, cluster_type == "IBL")
Seurat::Idents(subsobj) = subsobj$sample_type

table(subsobj$sample_type)

## 
##   HS   HD 
## 1284  226

We identify DE genes between HS and HD :

mark = Seurat::FindMarkers(subsobj, ident.1 = "HS", ident.2 = "HD")

mark = mark %>%
  dplyr::filter(p_val_adj < 0.05) %>%
  dplyr::arrange(-avg_logFC, pct.1 - pct.2)

list_results[["IBL_HS_vs_HD"]] = mark

dim(mark)

## [1] 107   5

head(mark, n = 20)

##                  p_val avg_logFC pct.1 pct.2    p_val_adj
## LGALS7    1.216786e-49 1.6340915 0.685 0.190 2.032154e-45
## S100A7    2.913857e-11 1.5790394 0.289 0.080 4.866433e-07
## MIF       1.195470e-44 1.2201587 0.706 0.261 1.996554e-40
## LGALS7B   1.099708e-09 1.1642854 0.769 0.717 1.836623e-05
## MTRNR2L12 9.867130e-12 1.0580493 0.393 0.181 1.647909e-07
## MTRNR2L8  3.047443e-11 0.9124631 0.350 0.155 5.089535e-07
## MT-CO1    4.350400e-12 0.9105978 0.449 0.226 7.265603e-08
## FOSB      7.804218e-09 0.8760606 0.298 0.133 1.303382e-04
## NDUFA4L2  1.283179e-20 0.8201600 0.743 0.434 2.143037e-16
## MT-CO2    1.272279e-08 0.8121340 0.433 0.265 2.124832e-04
## S100A9    1.727393e-13 0.7487053 0.315 0.075 2.884919e-09
## KRT15     2.873193e-08 0.6857105 0.290 0.115 4.798519e-04
## CA2       1.116064e-16 0.6640743 0.460 0.173 1.863938e-12
## ARF5      3.584966e-27 0.6221800 0.433 0.044 5.987251e-23
## CKB       4.051432e-08 0.5928841 0.282 0.115 6.766297e-04
## MT-ND4    2.524958e-07 0.5633098 0.410 0.248 4.216932e-03
## MT-ATP6   1.537778e-06 0.5539256 0.496 0.358 2.568244e-02
## CST6      1.707459e-07 0.5234863 0.421 0.221 2.851627e-03
## IGKC      5.879378e-07 0.5045779 0.101 0.000 9.819150e-03
## MCL1      1.177137e-07 0.5038666 0.379 0.208 1.965937e-03

ORS : HS vs HD

We subset the Seurat object and change cell identities to sample type.

subsobj = subset(sobj, cluster_type == "ORS")
Seurat::Idents(subsobj) = subsobj$sample_type

table(subsobj$sample_type)

## 
##   HS   HD 
## 1563  459

We identify DE genes between HS and HD :

mark = Seurat::FindMarkers(subsobj, ident.1 = "HS", ident.2 = "HD")

mark = mark %>%
  dplyr::filter(p_val_adj < 0.05) %>%
  dplyr::arrange(-avg_logFC, pct.1 - pct.2)

list_results[["ORS_HS_vs_HD"]] = mark

dim(mark)

## [1] 145   5

head(mark, n = 20)

##                   p_val avg_logFC pct.1 pct.2     p_val_adj
## S100A7     1.177379e-34 1.3981213 0.390 0.083  1.966341e-30
## S100A9     8.042694e-80 1.3882495 0.831 0.346  1.343210e-75
## S100A8     6.069068e-63 1.2071418 0.784 0.353  1.013595e-58
## RPS26     8.289007e-157 1.0156704 0.971 0.950 1.384347e-152
## LGALS7B    1.751974e-26 0.9144582 0.636 0.503  2.925971e-22
## CCL2       4.001592e-35 0.8469394 0.592 0.279  6.683058e-31
## LTF        1.036558e-17 0.7369496 0.183 0.022  1.731156e-13
## MTRNR2L8   2.332458e-41 0.7044528 0.909 0.904  3.895438e-37
## IFI27      4.002962e-19 0.6856359 0.194 0.024  6.685347e-15
## TIMP1      1.883812e-25 0.6120579 0.566 0.342  3.146155e-21
## IFITM3     2.955422e-49 0.5618141 0.889 0.819  4.935850e-45
## CD74       9.221913e-27 0.5594045 0.347 0.096  1.540152e-22
## HLA-C      1.036948e-24 0.5483711 0.731 0.547  1.731807e-20
## MIF        1.151305e-28 0.5462406 0.403 0.159  1.922795e-24
## AKR1B10    7.040926e-44 0.5197609 0.372 0.033  1.175905e-39
## MTRNR2L12  1.018482e-25 0.4788943 0.951 0.943  1.700967e-21
## MTRNR2L10  1.055574e-11 0.4753293 0.589 0.538  1.762914e-07
## CXCL14     2.007187e-14 0.4495196 0.882 0.802  3.352203e-10
## LGALS7     5.171905e-12 0.4358938 0.232 0.096  8.637599e-08
## IL1R2      9.715791e-13 0.4197815 0.226 0.085  1.622634e-08

Heatmap

We represent differentially expressed genes. First, we extract all DE genes :

features_oi = c(mark_to_label$gene_name,
                rownames(list_results$IBL_HS_vs_HD),
                rownames(list_results$ORS_HS_vs_HD)) %>%
  unique()

length(features_oi)

## [1] 259

Heatmap

We prepare the scaled expression matrix :

mat_expression = Seurat::GetAssayData(sobj, assay = "RNA", slot = "data")[features_oi, ]
mat_expression = Matrix::t(mat_expression)
mat_expression = dynutils::scale_quantile(mat_expression) # between 0 and 1
mat_expression = Matrix::t(mat_expression)
mat_expression = as.matrix(mat_expression) # not sparse

dim(mat_expression)

## [1]  259 3532

We prepare the heatmap annotation :

ha_top = ComplexHeatmap::HeatmapAnnotation(
  cluster_type = sobj$cluster_type,
  sample_type = sobj$sample_type,
  cluster = sobj$seurat_clusters,
  col = list(cluster_type = color_markers,
             sample_type = setNames(nm = c("HS", "HD"),
                                    c("#C55F40", "#2C78E6")),
             cluster = setNames(nm = levels(sobj$seurat_clusters),
                                aquarius::gg_color_hue(length(levels(sobj$seurat_clusters))))))

And the heatmap :

sobj$cell_group = paste0(sobj$cluster_type, sobj$sample_type) %>%
  as.factor()

ht = ComplexHeatmap::Heatmap(mat_expression,
                             col = aquarius::color_cnv,
                             # Annotation
                             top_annotation = ha_top,
                             # Grouping
                             column_order = sobj@meta.data %>%
                               dplyr::arrange(cluster_type, sample_type, seurat_clusters) %>%
                               rownames(),
                             column_split = sobj$cell_group,
                             column_gap = unit(c(0.01, 2, 0.01), "mm"),
                             column_title = NULL,
                             cluster_rows = TRUE,
                             cluster_columns = FALSE,
                             show_column_names = FALSE,
                             # Visual aspect
                             show_heatmap_legend = TRUE,
                             border = TRUE)

ComplexHeatmap::draw(ht,
                     merge_legend = TRUE,
                     heatmap_legend_side = "bottom",
                     annotation_legend_side = "bottom")

Save

We save the list of results :

saveRDS(list_results, file = paste0(out_dir, "/", save_name, "_list_results.rds"))

R Session

show

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/local/lib/R/lib/libRblas.so
## LAPACK: /usr/local/lib/R/lib/libRlapack.so
## 
## locale:
## [1] C
## 
## attached base packages:
##  [1] parallel  stats4    grid      stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] org.Mm.eg.db_3.10.0   AnnotationDbi_1.48.0  IRanges_2.20.2       
##  [4] S4Vectors_0.24.4      Biobase_2.46.0        BiocGenerics_0.32.0  
##  [7] ComplexHeatmap_2.14.0 ggplot2_3.3.5         patchwork_1.1.2      
## [10] dplyr_1.0.7          
## 
## loaded via a namespace (and not attached):
##   [1] softImpute_1.4              graphlayouts_0.7.0         
##   [3] pbapply_1.4-2               lattice_0.20-41            
##   [5] haven_2.3.1                 vctrs_0.3.8                
##   [7] usethis_2.0.1               dynwrap_1.2.1              
##   [9] blob_1.2.1                  survival_3.2-13            
##  [11] prodlim_2019.11.13          dynutils_1.0.5             
##  [13] later_1.3.0                 DBI_1.1.1                  
##  [15] R.utils_2.11.0              SingleCellExperiment_1.8.0 
##  [17] rappdirs_0.3.3              uwot_0.1.8                 
##  [19] dqrng_0.2.1                 jpeg_0.1-8.1               
##  [21] zlibbioc_1.32.0             pspline_1.0-18             
##  [23] pcaMethods_1.78.0           mvtnorm_1.1-1              
##  [25] htmlwidgets_1.5.4           GlobalOptions_0.1.2        
##  [27] future_1.22.1               UpSetR_1.4.0               
##  [29] laeken_0.5.2                leiden_0.3.3               
##  [31] clustree_0.4.3              scater_1.14.6              
##  [33] irlba_2.3.3                 DEoptimR_1.0-9             
##  [35] tidygraph_1.1.2             Rcpp_1.0.9                 
##  [37] readr_2.0.2                 KernSmooth_2.23-17         
##  [39] carrier_0.1.0               promises_1.1.0             
##  [41] gdata_2.18.0                DelayedArray_0.12.3        
##  [43] limma_3.42.2                graph_1.64.0               
##  [45] RcppParallel_5.1.4          Hmisc_4.4-0                
##  [47] fs_1.5.2                    RSpectra_0.16-0            
##  [49] fastmatch_1.1-0             ranger_0.12.1              
##  [51] digest_0.6.25               png_0.1-7                  
##  [53] sctransform_0.2.1           cowplot_1.0.0              
##  [55] DOSE_3.12.0                 here_1.0.1                 
##  [57] TInGa_0.0.0.9000            ggraph_2.0.3               
##  [59] pkgconfig_2.0.3             GO.db_3.10.0               
##  [61] DelayedMatrixStats_1.8.0    gower_0.2.1                
##  [63] ggbeeswarm_0.6.0            iterators_1.0.12           
##  [65] DropletUtils_1.6.1          reticulate_1.26            
##  [67] clusterProfiler_3.14.3      SummarizedExperiment_1.16.1
##  [69] circlize_0.4.15             beeswarm_0.4.0             
##  [71] GetoptLong_1.0.5            xfun_0.35                  
##  [73] bslib_0.3.1                 zoo_1.8-10                 
##  [75] tidyselect_1.1.0            reshape2_1.4.4             
##  [77] purrr_0.3.4                 ica_1.0-2                  
##  [79] pcaPP_1.9-73                viridisLite_0.3.0          
##  [81] rtracklayer_1.46.0          rlang_1.0.2                
##  [83] hexbin_1.28.1               jquerylib_0.1.4            
##  [85] dyneval_0.9.9               glue_1.4.2                 
##  [87] RColorBrewer_1.1-2          matrixStats_0.56.0         
##  [89] stringr_1.4.0               lava_1.6.7                 
##  [91] europepmc_0.3               DESeq2_1.26.0              
##  [93] recipes_0.1.17              labeling_0.3               
##  [95] httpuv_1.5.2                class_7.3-17               
##  [97] BiocNeighbors_1.4.2         DO.db_2.9                  
##  [99] annotate_1.64.0             jsonlite_1.7.2             
## [101] XVector_0.26.0              bit_4.0.4                  
## [103] mime_0.9                    aquarius_0.1.5             
## [105] Rsamtools_2.2.3             gridExtra_2.3              
## [107] gplots_3.0.3                stringi_1.4.6              
## [109] processx_3.5.2              gsl_2.1-6                  
## [111] bitops_1.0-6                cli_3.0.1                  
## [113] batchelor_1.2.4             RSQLite_2.2.0              
## [115] randomForest_4.6-14         tidyr_1.1.4                
## [117] data.table_1.14.2           rstudioapi_0.13            
## [119] GenomicAlignments_1.22.1    nlme_3.1-147               
## [121] qvalue_2.18.0               scran_1.14.6               
## [123] locfit_1.5-9.4              scDblFinder_1.1.8          
## [125] listenv_0.8.0               ggthemes_4.2.4             
## [127] gridGraphics_0.5-0          R.oo_1.24.0                
## [129] dbplyr_1.4.4                TTR_0.24.2                 
## [131] readxl_1.3.1                lifecycle_1.0.1            
## [133] timeDate_3043.102           ggpattern_0.3.1            
## [135] munsell_0.5.0               cellranger_1.1.0           
## [137] R.methodsS3_1.8.1           proxyC_0.1.5               
## [139] visNetwork_2.0.9            caTools_1.18.0             
## [141] codetools_0.2-16            GenomeInfoDb_1.22.1        
## [143] vipor_0.4.5                 lmtest_0.9-38              
## [145] msigdbr_7.5.1               htmlTable_1.13.3           
## [147] triebeard_0.3.0             lsei_1.2-0                 
## [149] xtable_1.8-4                ROCR_1.0-7                 
## [151] BiocManager_1.30.10         scatterplot3d_0.3-41       
## [153] abind_1.4-5                 farver_2.0.3               
## [155] parallelly_1.28.1           RANN_2.6.1                 
## [157] askpass_1.1                 GenomicRanges_1.38.0       
## [159] RcppAnnoy_0.0.16            tibble_3.1.5               
## [161] ggdendro_0.1-20             cluster_2.1.0              
## [163] future.apply_1.5.0          Seurat_3.1.5               
## [165] dendextend_1.15.1           Matrix_1.3-2               
## [167] ellipsis_0.3.2              prettyunits_1.1.1          
## [169] lubridate_1.7.9             ggridges_0.5.2             
## [171] igraph_1.2.5                RcppEigen_0.3.3.7.0        
## [173] fgsea_1.12.0                remotes_2.4.2              
## [175] scBFA_1.0.0                 destiny_3.0.1              
## [177] VIM_6.1.1                   testthat_3.1.0             
## [179] htmltools_0.5.2             BiocFileCache_1.10.2       
## [181] yaml_2.2.1                  utf8_1.1.4                 
## [183] plotly_4.9.2.1              XML_3.99-0.3               
## [185] ModelMetrics_1.2.2.2        e1071_1.7-3                
## [187] foreign_0.8-76              withr_2.5.0                
## [189] fitdistrplus_1.0-14         BiocParallel_1.20.1        
## [191] xgboost_1.4.1.1             bit64_4.0.5                
## [193] foreach_1.5.0               robustbase_0.93-9          
## [195] Biostrings_2.54.0           GOSemSim_2.13.1            
## [197] rsvd_1.0.3                  memoise_2.0.0              
## [199] evaluate_0.18               forcats_0.5.0              
## [201] rio_0.5.16                  geneplotter_1.64.0         
## [203] tzdb_0.1.2                  caret_6.0-86               
## [205] ps_1.6.0                    DiagrammeR_1.0.6.1         
## [207] curl_4.3                    fdrtool_1.2.15             
## [209] fansi_0.4.1                 highr_0.8                  
## [211] urltools_1.7.3              xts_0.12.1                 
## [213] GSEABase_1.48.0             acepack_1.4.1              
## [215] edgeR_3.28.1                checkmate_2.0.0            
## [217] scds_1.2.0                  cachem_1.0.6               
## [219] npsurv_0.4-0                babelgene_22.3             
## [221] rjson_0.2.20                openxlsx_4.1.5             
## [223] ggrepel_0.9.1               clue_0.3-60                
## [225] rprojroot_2.0.2             stabledist_0.7-1           
## [227] tools_3.6.3                 sass_0.4.0                 
## [229] nichenetr_1.1.1             magrittr_2.0.1             
## [231] RCurl_1.98-1.2              proxy_0.4-24               
## [233] car_3.0-11                  ape_5.3                    
## [235] ggplotify_0.0.5             xml2_1.3.2                 
## [237] httr_1.4.2                  assertthat_0.2.1           
## [239] rmarkdown_2.18              boot_1.3-25                
## [241] globals_0.14.0              R6_2.4.1                   
## [243] Rhdf5lib_1.8.0              nnet_7.3-14                
## [245] RcppHNSW_0.2.0              progress_1.2.2             
## [247] genefilter_1.68.0           statmod_1.4.34             
## [249] gtools_3.8.2                shape_1.4.6                
## [251] HDF5Array_1.14.4            BiocSingular_1.2.2         
## [253] rhdf5_2.30.1                splines_3.6.3              
## [255] AUCell_1.8.0                carData_3.0-4              
## [257] colorspace_1.4-1            generics_0.1.0             
## [259] base64enc_0.1-3             dynfeature_1.0.0           
## [261] smoother_1.1                gridtext_0.1.1             
## [263] pillar_1.6.3                tweenr_1.0.1               
## [265] sp_1.4-1                    ggplot.multistats_1.0.0    
## [267] rvcheck_0.1.8               GenomeInfoDbData_1.2.2     
## [269] plyr_1.8.6                  gtable_0.3.0               
## [271] zip_2.2.0                   knitr_1.41                 
## [273] latticeExtra_0.6-29         biomaRt_2.42.1             
## [275] fastmap_1.1.0               ADGofTest_0.3              
## [277] copula_1.0-0                doParallel_1.0.15          
## [279] vcd_1.4-8                   babelwhale_1.0.1           
## [281] openssl_1.4.1               scales_1.1.1               
## [283] backports_1.2.1             ipred_0.9-12               
## [285] enrichplot_1.6.1            hms_1.1.1                  
## [287] ggforce_0.3.1               Rtsne_0.15                 
## [289] shiny_1.7.1                 numDeriv_2016.8-1.1        
## [291] polyclip_1.10-0             lazyeval_0.2.2             
## [293] Formula_1.2-3               tsne_0.1-3                 
## [295] crayon_1.3.4                MASS_7.3-54                
## [297] pROC_1.16.2                 viridis_0.5.1              
## [299] dynparam_1.0.0              rpart_4.1-15               
## [301] zinbwave_1.8.0              compiler_3.6.3             
## [303] ggtext_0.1.0